library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 0.8.3 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(glmnet)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
##
## expand
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loaded glmnet 2.0-18
Note on how data are structured in each csv: First 30 are non drug Next 30 are drug
Alternate between left and right turning per row
Please be sure to put the data files under the /data folder
### Read in all the datafiles
r1<-(read.csv('data/r1.csv', header = FALSE))
r2<-(read.csv('data/r2.csv', header = FALSE))
r3<-(read.csv('data/r3.csv', header = FALSE))
r4<-(read.csv('data/r4.csv', header = FALSE))
r5<-(read.csv('data/r5.csv', header = FALSE))
r6<-(read.csv('data/r6.csv', header = FALSE))
r7<-(read.csv('data/r7.csv', header = FALSE))
r8<-(read.csv('data/r8.csv', header = FALSE))
## Create two vectors for outcome variables based on how we know the data are structured
drug_cond<-c(rep(0,30),rep(1,30))
dir<- c(rep(c(0,1),30))
y_colnames = c('drug_cond','dir')
### Function returns mean cross validated accuracy for glmnet models for each rat
### So 8 rats * 6 cv splits* 2 accuracy scores (drug condition/direction)
manual_cv<- function(df){
running_acc_dg = {}
running_acc_dir = {}
for (i in 0:5) {
df=cbind(df,drug_cond,dir)
s=(10*i)+1
e=10*(i+1)
test_set = c(s:e)
fit_df = df[!(rownames(df)%in%test_set),]
test_df = df[test_set,]
fit_dg = fit_df[,'drug_cond']
fit_dir = fit_df[,'dir']
test_dg = test_df[,'drug_cond']
test_dir = test_df[,'dir']
data_cols = as.matrix(fit_df[,!(colnames(fit_df)%in%y_colnames)])
test_cols = as.matrix(test_df[,!(colnames(test_df)%in%y_colnames)])
model_dg <- cv.glmnet(data_cols,fit_dg)
model_dir <- cv.glmnet(data_cols,fit_dir)
pred_dg = predict.cv.glmnet(model_dg,test_cols,s = "lambda.min")
pred_dir = predict.cv.glmnet(model_dir,test_cols,s = "lambda.min")
dg_acc= 1-(sum(abs(test_dg-pred_dg))/100)
dir_acc = 1 - (sum(abs(test_dir-pred_dir))/100)
running_acc_dg= append(running_acc_dg, dg_acc)
running_acc_dir= append(running_acc_dir, dir_acc)
plot(model_dg)
}
return_list = list(mean(running_acc_dg), mean(running_acc_dir))
return(return_list)
}
### running the function on all our
all_rats<- list(r1,r2,r3,r4,r5,r6,r7,r8)
j=0
for(this_rat in all_rats){
j=j+1
this_rat[is.na(this_rat)] <- 0
mean_acc<-manual_cv(this_rat)
print(paste0("mean drug prediction accuracy for rat number ",j," is ",mean_acc[1]))
print(paste0("mean direction prediction accuracy for rat number ",j," is ",mean_acc[2]))
}
## [1] "mean drug prediction accuracy for rat number 1 is 0.979039626237733"
## [1] "mean direction prediction accuracy for rat number 1 is 0.949600724794351"
## [1] "mean drug prediction accuracy for rat number 2 is 0.945991039580656"
## [1] "mean direction prediction accuracy for rat number 2 is 0.948509804974545"
## [1] "mean drug prediction accuracy for rat number 3 is 0.940122702928221"
## [1] "mean direction prediction accuracy for rat number 3 is 0.95149771951286"
## [1] "mean drug prediction accuracy for rat number 4 is 0.948218291046863"
## [1] "mean direction prediction accuracy for rat number 4 is 0.948193156404237"
## [1] "mean drug prediction accuracy for rat number 5 is 0.94067503014628"
## [1] "mean direction prediction accuracy for rat number 5 is 0.946832505854745"
## [1] "mean drug prediction accuracy for rat number 6 is 0.942569572329386"
## [1] "mean direction prediction accuracy for rat number 6 is 0.949935067589102"
## [1] "mean drug prediction accuracy for rat number 7 is 0.939184658223903"
## [1] "mean direction prediction accuracy for rat number 7 is 0.949994824428349"
## [1] "mean drug prediction accuracy for rat number 8 is 0.940983237596321"
## [1] "mean direction prediction accuracy for rat number 8 is 0.951103964746394"